#include "../common.h"

/*
 *  OpenCL version for two vector arrays A and B
 */
__kernel void corrcoefmatrixAB_opencl(__global float *out, __global float *inA, __global float *inB, const int n, const int m, const int tileSize) {
  __local float inAS[TILESIZE][TILESIZE];  // tileSize <= TILESIZE
  __local float inBS[TILESIZE][TILESIZE];

  int i,j,k;
  int tileBegin;
  float a_x,b_x;
  float sum_a, sum_b, sum_a2, sum_b2, sum_ab, corrcoef;

  i = get_global_id(1);
  j = get_global_id(0);

  sum_a = sum_a2 = sum_b = sum_b2 = sum_ab = 0;
  for (tileBegin=0; tileBegin < m; tileBegin += tileSize) {
     // Load data into shared memory
     inAS[get_local_id(1)][get_local_id(0)] = inA[get_global_id(1)*m+tileBegin+get_local_id(0)];
     inBS[get_local_id(1)][get_local_id(0)] = inB[(get_group_id(0)*get_local_size(0)+get_local_id(1))*m+tileBegin+get_local_id(0)];
     barrier(CLK_LOCAL_MEM_FENCE);

     // Do partial calculation on each tile
     for (k=0; k < tileSize; k++) {
        if (tileBegin + k < m) { // we should test if it is as good as padding
           a_x = inAS[get_local_id(1)][k];
           b_x = inBS[get_local_id(0)][k];
           sum_a += a_x;
           sum_a2 += a_x*a_x;
           sum_b += b_x;
           sum_b2 += b_x*b_x;
           sum_ab += a_x*b_x;
        }
     }
     barrier(CLK_LOCAL_MEM_FENCE);
  }
  corrcoef = (m*sum_ab - sum_a*sum_b)/sqrt((m*sum_a2-sum_a*sum_a)*(m*sum_b2-sum_b*sum_b));
  out[i*n+j] = corrcoef;
}

